## title:   Narcissism in Political Participation
## authors: Z. Fazekas & P.K. Hatemi
## study:   all
## goal:    run regression models, format output data for plots
##          export regression output in table format (for appendix)
## ----
## 0. packages and functions --
library("dplyr")
library("broom")
library("reshape2")
library("ggplot2")
library("ggthemes")
library("psych")
library("polycor")
library("texreg")

## 1. load data -- 
source("helper.R")
dk11 <- read.csv("./data/dk11-reg.csv")
us13 <- read.csv("./data/us13-reg.csv")
us15 <- read.csv("./data/us15-reg.csv")

## housekeeping to facilitate running all models in one shot
## see 1-descriptives_all.R as well
dk11$not_caucasian <- 0
dk11$exhib_sd      <- 0  
dk11$exhib         <- NA 
us15$eff_full      <- 0
names(dk11)[1]     <- "uid"
names(us13)[1]     <- "uid"
names(us15)[1]     <- "uid"

us13$sample_num    <- 1
us13$sample_num[us13$sample == "Wave2"]    <- 2
us13$sample_num[us13$sample == "Oversample"]    <- 3
dk11$sample_num    <- 1
us15$sample_num    <- 1

us13$weight <- 1
dk11$weight <- 1

## 2. creating long format + merged data for regression models

dk11$interest <- round(dk11$interest*3, 0)
us13$interest <- round(us13$interest*3, 0)
us15$interest <- round(us15$interest*3, 0)

dk11$interest_sd <- (dk11$interest - mean(dk11$interest, na.rm = TRUE))/
  2*sd(dk11$interest, na.rm = TRUE)
us13$interest_sd <- (us13$interest - mean(us13$interest, na.rm = TRUE))/
  2*sd(us13$interest, na.rm = TRUE)
us15$interest_sd <- (us15$interest - mean(us15$interest, na.rm = TRUE))/
  2*sd(us15$interest, na.rm = TRUE)

outcomes <- c("part_full")
id_vars  <- c("uid", "sample_num", "study", "weight", 
              "female", "edu_cat", "age_sd", "not_caucasian",
              "interest_sd",
              "npi_sd",
              "authority_sd", "exhib_sd", "superior_sd", "exploit_sd",
              "entitle_sd", "vanity_sd", "sufficient_sd",
              "npi25_sd", "leadauth_sd",
              "entexp_sd", "grandexhib_sd")

all_vars <- c(id_vars, outcomes)

all_reg <- rbind(dk11[, all_vars],
                 us13[, all_vars],
                 us15[, all_vars])

all_reg_long <- melt(all_reg, id = id_vars)
to_rename <- (dim(all_reg_long)[2] - 1):(dim(all_reg_long)[2]) ## last 2 vars
names(all_reg_long)[to_rename] <- c("outcome", "y")

## models for continuous outcome
npi.mods <- all_reg_long %>%
            group_by(study, outcome) %>%
            do(npi_model = lm(y ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian, 
                                  weights = weight,
                                  data = .))

npi.mods.all <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(npi_model = lm(y ~ npi_sd + 
                      female + age_sd + edu_cat +
                      not_caucasian + interest_sd, 
                    weights = weight,
                    data = .))

## run the 3factor models
npi25.mods <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(npi_model = lm(y ~ npi25_sd + 
                      female + age_sd + edu_cat +
                      not_caucasian, 
                    weights = weight,
                    data = .))
npi25.mods.all <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(npi_model = lm(y ~ npi25_sd + 
                      female + age_sd + edu_cat +
                      not_caucasian + interest_sd, 
                    weights = weight,
                    data = .))

## all subfacets as predictors
sub.mods <- all_reg_long %>%
            group_by(study, outcome) %>%
            do(facet_model = lm(y ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian, 
                                  weights = weight,
                                  data = .))
sub.mods.all <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(facet_model = lm(y ~ authority_sd +
                        superior_sd +
                        exploit_sd +
                        entitle_sd +
                        sufficient_sd +
                        vanity_sd +
                        exhib_sd + 
                        female + age_sd + edu_cat +
                        not_caucasian + interest_sd, 
                      weights = weight,
                      data = .))

sub25.mods <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(facet_model = lm(y ~ leadauth_sd + entexp_sd + grandexhib_sd +
                        female + age_sd + edu_cat +
                        not_caucasian, 
                      weights = weight,
                      data = .))
sub25.mods.all <- all_reg_long %>%
  group_by(study, outcome) %>%
  do(facet_model = lm(y ~ leadauth_sd + entexp_sd + grandexhib_sd +
                        female + age_sd + edu_cat +
                        not_caucasian + interest_sd, 
                      weights = weight,
                      data = .))

## check and extract
npi_coefs <- npi.mods %>% tidy(npi_model) %>% data.frame()
npi25_coefs <- npi25.mods %>% tidy(npi_model) %>% data.frame()

npi_all_coefs <- npi.mods.all %>% tidy(npi_model) %>% data.frame()
npi25_all_coefs <- npi25.mods.all %>% tidy(npi_model) %>% data.frame()


npi_stats <- npi.mods %>% glance(npi_model) %>% data.frame()
npi25_stats <- npi25.mods %>% glance(npi_model) %>% data.frame()

npi_all_stats <- npi.mods.all %>% glance(npi_model) %>% data.frame()
npi25_all_stats <- npi25.mods.all %>% glance(npi_model) %>% data.frame()

## drop those that were constants (not included in the some studies)
npi_coefs <- dplyr::filter(npi_coefs, p.value   != "NaN")
npi25_coefs <- dplyr::filter(npi25_coefs, p.value   != "NaN")
npi_all_coefs <- dplyr::filter(npi_all_coefs, p.value   != "NaN")
npi25_all_coefs <- dplyr::filter(npi25_all_coefs, p.value   != "NaN")

npi_stats <- dplyr::filter(npi_stats, r.squared != "NaN")
npi_all_stats <- dplyr::filter(npi_all_stats, r.squared != "NaN")

npi25_stats <- dplyr::filter(npi25_stats, r.squared != "NaN")
npi25_all_stats <- dplyr::filter(npi25_all_stats, r.squared != "NaN")

npi_coefs$model_id <- paste(npi_coefs$study,
                            npi_coefs$outcome,
                            sep = "-")
npi_stats$model_id <- paste(npi_stats$study,
                            npi_stats$outcome,
                            sep = "-")
npi_all_coefs$model_id <- paste(npi_all_coefs$study,
                                npi_all_coefs$outcome,
                                sep = "-")
npi_all_stats$model_id <- paste(npi_all_stats$study,
                                npi_all_stats$outcome,
                                sep = "-")

## better labeling
npi_stats$outcome_str <- "Participation"

## adjusted R2
npi_stats$adj_r2 <- paste0("adj. R2 = ",
                           round(npi_stats$adj.r.squared, 3))
## sample size (for models)
npi_stats$n <- paste0("N = ", npi_stats$df + npi_stats$df.residual)
npi_coefs <- merge(npi_coefs, 
                   npi_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                   by = "model_id")
npi_all_stats$outcome_str <- "Participation"

## adjusted R2
npi_all_stats$adj_r2 <- paste0("adj. R2 = ",
                               round(npi_all_stats$adj.r.squared, 3))
## sample size (for models)
npi_all_stats$n <- paste0("N = ", npi_all_stats$df + npi_all_stats$df.residual)
npi_all_coefs <- merge(npi_all_coefs, 
                       npi_all_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                       by = "model_id")

## idem for models using the subfacets
sub_coefs <- sub.mods %>% tidy(facet_model) %>% data.frame()
sub_stats <- sub.mods %>% glance(facet_model) %>% data.frame()

sub_coefs <- dplyr::filter(sub_coefs, p.value   != "NaN")
sub_stats <- dplyr::filter(sub_stats, r.squared != "NaN")

sub_coefs$model_id <- paste(sub_coefs$study,
                            sub_coefs$outcome,
                            sep = "-")
sub_stats$model_id <- paste(sub_stats$study,
                            sub_stats$outcome,
                            sep = "-")
sub_stats$outcome_str <- "Participation"

sub_stats$adj_r2 <- paste0("adj. R2 = ",
                           round(sub_stats$adj.r.squared, 3))
sub_stats$n <- paste0("N = ", sub_stats$df + sub_stats$df.residual)
sub_coefs <- merge(sub_coefs, 
                   sub_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                   by = "model_id")

## add model name (full smmed scale or subfacets to each res. data)
npi_coefs$model <- "summed_scale"
sub_coefs$model <- "sub_facets"

## idem for models using the subfacets
sub_all_coefs <- sub.mods.all %>% tidy(facet_model) %>% data.frame()
sub_all_stats <- sub.mods.all %>% glance(facet_model) %>% data.frame()

sub_all_coefs <- dplyr::filter(sub_all_coefs, p.value   != "NaN")
sub_all_stats <- dplyr::filter(sub_all_stats, r.squared != "NaN")

sub_all_coefs$model_id <- paste(sub_all_coefs$study,
                                sub_all_coefs$outcome,
                                sep = "-")
sub_all_stats$model_id <- paste(sub_all_stats$study,
                                sub_all_stats$outcome,
                                sep = "-")
sub_all_stats$outcome_str <- "Participation"

sub_all_stats$adj_r2 <- paste0("adj. R2 = ",
                               round(sub_all_stats$adj.r.squared, 3))
sub_all_stats$n <- paste0("N = ", sub_all_stats$df + sub_all_stats$df.residual)
sub_all_coefs <- merge(sub_all_coefs, 
                       sub_all_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                       by = "model_id")

## add model name (full smmed scale or subfacets to each res. data)
npi_all_coefs$model <- "summed_scale"
sub_all_coefs$model <- "sub_facets"

npi_all_coefs$spec <- "all"
sub_all_coefs$spec <- "all"

npi_coefs$spec <- "socdem"
sub_coefs$spec <- "socdem"

model_results <- rbind(npi_coefs, sub_coefs,
                       npi_all_coefs, sub_all_coefs)

## labeling terms
model_results$term_str <- "Narcissism"
model_results$term_str[model_results$term == "(Intercept)"] <- 
    "Intercept"
model_results$term_str[model_results$term == "female"] <- 
    "Gender (Female = 1)"
model_results$term_str[model_results$term == "age_sd"] <- 
    "Age (2SD stand)"
model_results$term_str[model_results$term == "disc_sd"] <- 
  "Discussion (2SD stand)"
model_results$term_str[model_results$term == "interest_sd"] <- 
  "Interest (2SD stand)"
model_results$term_str[model_results$term == "edu_cat"] <- 
    "Higher education"
model_results$term_str[model_results$term == "not_caucasian"] <- 
    "Ethnicity/Race (Not white = 1)"
model_results$term_str[model_results$term == "not_caucasian"] <- 
    "Ethnicity/Race (Not white = 1)"
model_results$term_str[model_results$term == "authority_sd"] <- 
    "Authority"
model_results$term_str[model_results$term == "entitle_sd"] <- 
    "Entitlement"
model_results$term_str[model_results$term == "exhib_sd"] <- 
    "Exhibitionism"
model_results$term_str[model_results$term == "exploit_sd"] <- 
    "Exploitativeness"
model_results$term_str[model_results$term == "superior_sd"] <- 
    "Superiority"
model_results$term_str[model_results$term == "sufficient_sd"] <- 
    "Self-sufficiency"
model_results$term_str[model_results$term == "vanity_sd"] <- 
    "Vanity"


npi25_coefs$model_id <- paste(npi25_coefs$study,
                              npi25_coefs$outcome,
                              sep = "-")
npi25_stats$model_id <- paste(npi25_stats$study,
                              npi25_stats$outcome,
                              sep = "-")
npi25_all_coefs$model_id <- paste(npi25_all_coefs$study,
                                  npi25_all_coefs$outcome,
                                  sep = "-")
npi25_all_stats$model_id <- paste(npi25_all_stats$study,
                                  npi25_all_stats$outcome,
                                  sep = "-")

## better labeling
npi25_stats$outcome_str <- "Participation"

## adjusted R2
npi25_stats$adj_r2 <- paste0("adj. R2 = ",
                             round(npi25_stats$adj.r.squared, 3))
## sample size (for models)
npi25_stats$n <- paste0("N = ", npi25_stats$df + npi25_stats$df.residual)
npi25_coefs <- merge(npi25_coefs, 
                     npi25_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                     by = "model_id")
npi25_all_stats$outcome_str <- "Participation"

## adjusted R2
npi25_all_stats$adj_r2 <- paste0("adj. R2 = ",
                                 round(npi25_all_stats$adj.r.squared, 3))
## sample size (for models)
npi25_all_stats$n <- paste0("N = ", npi25_all_stats$df + npi25_all_stats$df.residual)
npi25_all_coefs <- merge(npi25_all_coefs, 
                         npi25_all_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                         by = "model_id")

## idem for models using the subfacets
sub25_coefs <- sub25.mods %>% tidy(facet_model) %>% data.frame()
sub25_stats <- sub25.mods %>% glance(facet_model) %>% data.frame()

sub25_coefs <- dplyr::filter(sub25_coefs, p.value   != "NaN")
sub25_stats <- dplyr::filter(sub25_stats, r.squared != "NaN")

sub25_coefs$model_id <- paste(sub25_coefs$study,
                              sub25_coefs$outcome,
                              sep = "-")
sub25_stats$model_id <- paste(sub25_stats$study,
                              sub25_stats$outcome,
                              sep = "-")
sub25_stats$outcome_str <- "Participation"

sub25_stats$adj_r2 <- paste0("adj. R2 = ",
                             round(sub25_stats$adj.r.squared, 3))
sub25_stats$n <- paste0("N = ", sub25_stats$df + sub25_stats$df.residual)
sub25_coefs <- merge(sub25_coefs, 
                     sub25_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                     by = "model_id")

## add model name (full smmed scale or subfacets to each res. data)
npi25_coefs$model <- "summed_scale"
sub25_coefs$model <- "sub_facets"

## idem for models using the subfacets
sub25_all_coefs <- sub25.mods.all %>% tidy(facet_model) %>% data.frame()
sub25_all_stats <- sub25.mods.all %>% glance(facet_model) %>% data.frame()

sub25_all_coefs <- dplyr::filter(sub25_all_coefs, p.value   != "NaN")
sub25_all_stats <- dplyr::filter(sub25_all_stats, r.squared != "NaN")

sub25_all_coefs$model_id <- paste(sub25_all_coefs$study,
                                  sub25_all_coefs$outcome,
                                  sep = "-")
sub25_all_stats$model_id <- paste(sub25_all_stats$study,
                                  sub25_all_stats$outcome,
                                  sep = "-")
sub25_all_stats$outcome_str <- "Participation"

sub25_all_stats$adj_r2 <- paste0("adj. R2 = ",
                                 round(sub25_all_stats$adj.r.squared, 3))
sub25_all_stats$n <- paste0("N = ", sub25_all_stats$df + sub25_all_stats$df.residual)
sub25_all_coefs <- merge(sub25_all_coefs, 
                         sub25_all_stats[, c("model_id", "outcome_str", "adj_r2", "n")],
                         by = "model_id")

## add model name (full smmed scale or subfacets to each res. data)
npi25_all_coefs$model <- "summed_scale"
sub25_all_coefs$model <- "sub_facets"

npi25_all_coefs$spec <- "all"
sub25_all_coefs$spec <- "all"

npi25_coefs$spec <- "socdem"
sub25_coefs$spec <- "socdem"

model25_results <- rbind(npi25_coefs, sub25_coefs,
                         npi25_all_coefs, sub25_all_coefs)

model25_results$term_str <- "Narcissism (25)"
model25_results$term_str[model25_results$term == "(Intercept)"] <- 
  "Intercept"
model25_results$term_str[model25_results$term == "female"] <- 
  "Gender (Female = 1)"
model25_results$term_str[model25_results$term == "age_sd"] <- 
  "Age (2SD stand)"
model25_results$term_str[model25_results$term == "disc_sd"] <- 
  "Discussion (2SD stand)"
model25_results$term_str[model25_results$term == "interest_sd"] <- 
  "Interest (2SD stand)"
model25_results$term_str[model25_results$term == "edu_cat"] <- 
  "Higher education"
model25_results$term_str[model25_results$term == "not_caucasian"] <- 
  "Ethnicity/Race (Not white = 1)"
model25_results$term_str[model25_results$term == "not_caucasian"] <- 
  "Ethnicity/Race (Not white = 1)"
model25_results$term_str[model25_results$term == "leadauth_sd"] <- 
  "Leadership/Authority"
model25_results$term_str[model25_results$term == "entexp_sd"] <- 
  "Entitlement/Exploitativeness"
model25_results$term_str[model25_results$term == "grandexhib_sd"] <- 
  "Grandiose/Exhibitionism"

model_results$narc_op <- "Full"
model25_results$narc_op <- "Three"


## 3. export datasets --
all_mods <- rbind(model_results, model25_results)
# write.csv(all_mods, file = "./data/model-coef-all.csv",
#           row.names = FALSE)

## 4. turnout models
m_turnout_dk11_npi   <- glm(turnout ~ npi_sd + 
                              female + age_sd + edu_cat,
                            data = dk11,
                            family = binomial(link = "logit"))
summary(m_turnout_dk11_npi)
m_turnout_dk11_facet <- glm(turnout ~ authority_sd +
                              superior_sd +
                              exploit_sd +
                              entitle_sd +
                              sufficient_sd +
                              vanity_sd +
                              exhib_sd + 
                              female + age_sd + edu_cat,
                            data = dk11,
                            family = binomial(link = "logit"))
summary(m_turnout_dk11_facet)
m_turnout_dk11_all_npi   <- glm(turnout ~ npi_sd + 
                                  female + age_sd + edu_cat + interest_sd,
                                data = dk11,
                                family = binomial(link = "logit"))
summary(m_turnout_dk11_all_npi)
m_turnout_dk11_all_facet <- glm(turnout ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat + interest_sd,
                                data = dk11,
                                family = binomial(link = "logit"))
summary(m_turnout_dk11_all_facet)


m_turnout_us13_npi   <- glm(turnout_2012 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us13,
                                  family = binomial(link = "logit"))
summary(m_turnout_us13_npi)
m_turnout_us13_facet <- glm(turnout_2012 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us13,
                                  family = binomial(link = "logit"))
summary(m_turnout_us13_facet)
m_turnout_us13_all_npi   <- glm(turnout_2012 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us13,
                                family = binomial(link = "logit"))
summary(m_turnout_us13_all_npi)
m_turnout_us13_all_facet <- glm(turnout_2012 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us13,
                                family = binomial(link = "logit"))
summary(m_turnout_us13_all_facet)

m_turnout_us15_npi   <- glm(turnout_2012 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))
summary(m_turnout_us15_npi)
m_turnout_us15_facet <- glm(turnout_2012 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))

m_midterm_us15_npi   <- glm(turnout_2014 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))
m_midterm_us15_facet <- glm(turnout_2014 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))
summary(m_midterm_us15_facet)

m_turnout_us15_all_npi   <- glm(turnout_2012 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))
m_turnout_us15_all_facet <- glm(turnout_2012 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))

m_midterm_us15_all_npi  <- glm(turnout_2014 ~ npi_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))
m_midterm_us15_all_facet <- glm(turnout_2014 ~ authority_sd +
                                  superior_sd +
                                  exploit_sd +
                                  entitle_sd +
                                  sufficient_sd +
                                  vanity_sd +
                                  exhib_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))


turnout_res <- rbind(
cbind(filter(tidy(m_turnout_dk11_npi), term == "npi_sd"), 
        outcome = "turnout", study = "DK11"),
cbind(filter(tidy(m_turnout_us13_npi), term == "npi_sd"), 
            outcome = "turnout_2012", study = "US13"),
cbind(filter(tidy(m_turnout_us15_npi), term == "npi_sd"), 
            outcome = "turnout_2012", study = "US15"),
cbind(filter(tidy(m_midterm_us15_npi), term == "npi_sd"), 
            outcome = "turnout_2014", study = "US15"),
cbind(tidy(m_turnout_dk11_facet)[2:7,],
      outcome = "turnout", study = "DK11"),
cbind(tidy(m_turnout_us13_facet)[2:8,],
            outcome = "turnout_2012", study = "US13"),
cbind(tidy(m_midterm_us15_facet)[2:8,],
            outcome = "turnout_2014", study = "US15"),
cbind(tidy(m_turnout_us15_facet)[2:8,],
            outcome = "turnout_2012", study = "US15")
)

turnout_res_all <- rbind(
  cbind(filter(tidy(m_turnout_dk11_all_npi), term == "npi_sd"), 
        outcome = "turnout", study = "DK11"),
  cbind(filter(tidy(m_turnout_us13_all_npi), term == "npi_sd"), 
        outcome = "turnout_2012", study = "US13"),
  cbind(filter(tidy(m_turnout_us15_all_npi), term == "npi_sd"), 
        outcome = "turnout_2012", study = "US15"),
  cbind(filter(tidy(m_midterm_us15_all_npi), term == "npi_sd"), 
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout_dk11_all_facet)[2:7,],
        outcome = "turnout", study = "DK11"),
  cbind(tidy(m_turnout_us13_all_facet)[2:8,],
        outcome = "turnout_2012", study = "US13"),
  cbind(tidy(m_midterm_us15_all_facet)[2:8,],
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout_us15_all_facet)[2:8,],
        outcome = "turnout_2012", study = "US15")
)

turnout_res$term_str <- "Narcissism"
turnout_res$term_str[turnout_res$term == "authority_sd"] <- 
    "Authority"
turnout_res$term_str[turnout_res$term == "entitle_sd"] <- 
    "Entitlement"
turnout_res$term_str[turnout_res$term == "exhib_sd"] <- 
    "Exhibitionism"
turnout_res$term_str[turnout_res$term == "exploit_sd"] <- 
    "Exploitativeness"
turnout_res$term_str[turnout_res$term == "superior_sd"] <- 
    "Superiority"
turnout_res$term_str[turnout_res$term == "sufficient_sd"] <- 
    "Self-sufficiency"
turnout_res$term_str[turnout_res$term == "vanity_sd"] <- 
    "Vanity"

turnout_res$outcome_str[turnout_res$outcome == "turnout_2012"] <- 
    "Turnout (2012 presidential)"
turnout_res$outcome_str[turnout_res$outcome == "turnout_2014"] <- 
    "Turnout (2014 midterm)"
turnout_res$outcome_str[turnout_res$outcome == "turnout"] <- 
  "Turnout (parliament)"

turnout_res_all$term_str <- "Narcissism"
turnout_res_all$term_str[turnout_res_all$term == "authority_sd"] <- 
  "Authority"
turnout_res_all$term_str[turnout_res_all$term == "entitle_sd"] <- 
  "Entitlement"
turnout_res_all$term_str[turnout_res_all$term == "exhib_sd"] <- 
  "Exhibitionism"
turnout_res_all$term_str[turnout_res_all$term == "exploit_sd"] <- 
  "Exploitativeness"
turnout_res_all$term_str[turnout_res_all$term == "superior_sd"] <- 
  "Superiority"
turnout_res_all$term_str[turnout_res_all$term == "sufficient_sd"] <- 
  "Self-sufficiency"
turnout_res_all$term_str[turnout_res_all$term == "vanity_sd"] <- 
  "Vanity"

turnout_res_all$outcome_str[turnout_res_all$outcome == "turnout_2012"] <- 
  "Turnout (2012 presidential)"
turnout_res_all$outcome_str[turnout_res_all$outcome == "turnout_2014"] <- 
  "Turnout (2014 midterm)"
turnout_res_all$outcome_str[turnout_res_all$outcome == "turnout"] <- 
  "Turnout (parliament)"

turnout_res$spec <- "socdem"
turnout_res_all$spec <- "all"

turnout_all <- rbind(turnout_res, turnout_res_all)

# 3-factor models
m_turnout25_dk11_npi   <- glm(turnout ~ npi25_sd + 
                                female + age_sd + edu_cat,
                              data = dk11,
                              family = binomial(link = "logit"))
summary(m_turnout25_dk11_npi)
m_turnout25_dk11_facet <- glm(turnout ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                female + age_sd + edu_cat,
                              data = dk11,
                              family = binomial(link = "logit"))
summary(m_turnout25_dk11_facet)
m_turnout25_dk11_all_npi   <- glm(turnout ~ npi25_sd + 
                                    female + age_sd + edu_cat + interest_sd,
                                  data = dk11,
                                  family = binomial(link = "logit"))
summary(m_turnout25_dk11_all_npi)
m_turnout25_dk11_all_facet <- glm(turnout ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                    female + age_sd + edu_cat + interest_sd,
                                  data = dk11,
                                  family = binomial(link = "logit"))
summary(m_turnout25_dk11_all_facet)


m_turnout25_us13_npi   <- glm(turnout_2012 ~ npi25_sd + 
                                female + age_sd + edu_cat +
                                not_caucasian,
                              data = us13,
                              family = binomial(link = "logit"))
summary(m_turnout25_us13_npi)
m_turnout25_us13_facet <- glm(turnout_2012 ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                female + age_sd + edu_cat +
                                not_caucasian,
                              data = us13,
                              family = binomial(link = "logit"))
summary(m_turnout25_us13_facet)
m_turnout25_us13_all_npi   <- glm(turnout_2012 ~ npi25_sd + 
                                    female + age_sd + edu_cat +
                                    not_caucasian + interest_sd,
                                  data = us13,
                                  family = binomial(link = "logit"))
summary(m_turnout25_us13_all_npi)
m_turnout25_us13_all_facet <- glm(turnout_2012 ~ leadauth_sd + entexp_sd + grandexhib_sd + 
                                    female + age_sd + edu_cat +
                                    not_caucasian + interest_sd,
                                  data = us13,
                                  family = binomial(link = "logit"))
summary(m_turnout25_us13_all_facet)


m_turnout25_us15_npi   <- glm(turnout_2012 ~ npi25_sd + 
                                female + age_sd + edu_cat +
                                not_caucasian,
                              data = us15,
                              # weights = weight,
                              family = binomial(link = "logit"))
summary(m_turnout25_us15_npi)
m_turnout25_us15_facet <- glm(turnout_2012 ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                female + age_sd + edu_cat +
                                not_caucasian,
                              data = us15,
                              # weights = weight,
                              family = binomial(link = "logit"))

m_midterm25_us15_npi   <- glm(turnout_2014 ~ npi25_sd + 
                              female + age_sd + edu_cat +
                              not_caucasian,
                            data = us15,
                            # weights = weight,
                            family = binomial(link = "logit"))
m_midterm25_us15_facet <- glm(turnout_2014 ~ leadauth_sd + entexp_sd + grandexhib_sd +
                              female + age_sd + edu_cat +
                              not_caucasian,
                            data = us15,
                            # weights = weight,
                            family = binomial(link = "logit"))

m_turnout25_us15_all_npi   <- glm(turnout_2012 ~ npi25_sd + 
                                    female + age_sd + edu_cat +
                                    not_caucasian + interest_sd,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))
m_turnout25_us15_all_facet <- glm(turnout_2012 ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                    female + age_sd + edu_cat +
                                    not_caucasian + interest_sd,
                                  data = us15,
                                  # weights = weight,
                                  family = binomial(link = "logit"))

m_midterm25_us15_all_npi   <- glm(turnout_2014 ~ npi25_sd + 
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))
summary(m_midterm25_us15_all_npi)
m_midterm25_us15_all_facet <- glm(turnout_2014 ~ leadauth_sd + entexp_sd + grandexhib_sd +
                                  female + age_sd + edu_cat +
                                  not_caucasian + interest_sd,
                                data = us15,
                                # weights = weight,
                                family = binomial(link = "logit"))


turnout25_res <- rbind(
  cbind(filter(tidy(m_turnout25_dk11_npi), term == "npi25_sd"), 
        outcome = "turnout", study = "DK11"),
  cbind(filter(tidy(m_turnout25_us13_npi), term == "npi25_sd"), 
        outcome = "turnout_2012", study = "US13"),
  cbind(filter(tidy(m_turnout25_us15_npi), term == "npi25_sd"), 
        outcome = "turnout_2012", study = "US15"),
  cbind(filter(tidy(m_midterm25_us15_npi), term == "npi25_sd"), 
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout25_dk11_facet)[2:4,],
        outcome = "turnout", study = "DK11"),
  cbind(tidy(m_turnout25_us13_facet)[2:4,],
        outcome = "turnout_2012", study = "US13"),
  cbind(tidy(m_midterm25_us15_facet)[2:4,],
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout25_us15_facet)[2:4,],
        outcome = "turnout_2012", study = "US15")
)


turnout25_res_all <- rbind(
  cbind(filter(tidy(m_turnout25_dk11_all_npi), term == "npi25_sd"), 
        outcome = "turnout", study = "DK11"),
  cbind(filter(tidy(m_turnout25_us13_all_npi), term == "npi25_sd"), 
        outcome = "turnout_2012", study = "US13"),
  cbind(filter(tidy(m_turnout25_us15_all_npi), term == "npi25_sd"), 
        outcome = "turnout_2012", study = "US15"),
  cbind(filter(tidy(m_midterm25_us15_all_npi), term == "npi25_sd"), 
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout25_dk11_all_facet)[2:4,],
        outcome = "turnout", study = "DK11"),
  cbind(tidy(m_turnout25_us13_all_facet)[2:4,],
        outcome = "turnout_2012", study = "US13"),
  cbind(tidy(m_midterm25_us15_all_facet)[2:4,],
        outcome = "turnout_2014", study = "US15"),
  cbind(tidy(m_turnout25_us15_all_facet)[2:4,],
        outcome = "turnout_2012", study = "US15")
)

turnout25_res$term_str <- "Narcissism (25)"
turnout25_res$term_str[turnout25_res$term == "leadauth_sd"] <- 
  "Leadership/Authority"
turnout25_res$term_str[turnout25_res$term == "entexp_sd"] <- 
  "Entitlement/Exploitativeness"
turnout25_res$term_str[turnout25_res$term == "grandexhib_sd"] <- 
  "Grandiose/Exhibitionism"


turnout25_res$outcome_str[turnout25_res$outcome == "turnout_2012"] <- 
  "Turnout (2012 presidential)"
turnout25_res$outcome_str[turnout25_res$outcome == "turnout_2014"] <- 
  "Turnout (2014 midterm)"
turnout25_res$outcome_str[turnout25_res$outcome == "turnout"] <- 
  "Turnout (parliament)"

turnout25_res_all$term_str <- "Narcissism (25)"
turnout25_res_all$term_str[turnout25_res_all$term == "leadauth_sd"] <- 
  "Leadership/Authority"
turnout25_res_all$term_str[turnout25_res_all$term == "entexp_sd"] <- 
  "Entitlement/Exploitativeness"
turnout25_res_all$term_str[turnout25_res_all$term == "grandexhib_sd"] <- 
  "Grandiose/Exhibitionism"

turnout25_res_all$outcome_str[turnout25_res_all$outcome == "turnout_2012"] <- 
  "Turnout (2012 presidential)"
turnout25_res_all$outcome_str[turnout25_res_all$outcome == "turnout_2014"] <- 
  "Turnout (2014 midterm)"
turnout25_res_all$outcome_str[turnout25_res_all$outcome == "turnout"] <- 
  "Turnout (parliament)"

turnout25_res$spec <- "socdem"
turnout25_res_all$spec <- "all"

turnout25_all <- rbind(turnout25_res, turnout25_res_all)

turnout_all$narc_op <- "Full"
turnout25_all$narc_op <- "Three"

## 5. export turnout dataset --
turnout_full <- rbind(turnout_all, turnout25_all)
# write.csv(turnout_full, file = "./data/turnout-coef-all.csv",
#           row.names = FALSE)

# Table exports
screenreg(list(npi.mods$npi_model[[1]],
             m_turnout_dk11_npi,
             npi.mods$npi_model[[2]],
             m_turnout_us13_npi,
             npi.mods$npi_model[[3]],
             m_turnout_us15_npi,
             m_midterm_us15_npi),
          custom.model.names = c("DK11 Participation",
                                 "DK11 Turnout",
                                 "US13 Participation",
                                 "US13 Turnout",
                                 "US15 Participation",
                                 "US15 Turnout",
                                 "US15 Midterm"),
          custom.coef.names = c("Intercept", "Narcissism (full)",
                                "Female", "Age (2SD)", "Higher education", "Not Caucasian"),
          caption = "(Appendix) Narcissism, political participation, and turnout",
          caption.above = TRUE)


screenreg(list(npi25.mods$npi_model[[1]],
               m_turnout25_dk11_npi,
               npi25.mods$npi_model[[2]],
               m_turnout25_us13_npi,
               npi25.mods$npi_model[[3]],
               m_turnout25_us15_npi,
               m_midterm25_us15_npi),
          custom.model.names = c("DK11 Participation",
                                 "DK11 Turnout",
                                 "US13 Participation",
                                 "US13 Turnout",
                                 "US15 Participation",
                                 "US15 Turnout",
                                 "US15 Midterm"),
          custom.coef.names = c("Intercept", "Narcissism (subset)",
                                "Female", "Age (2SD)", "Higher education", "Not Caucasian"),
          caption = "(Appendix) Narcissism, political participation, and turnout (25 NPI items)",
          caption.above = TRUE)

# Same with interest
screenreg(list(npi.mods.all$npi_model[[1]],
             m_turnout_dk11_all_npi,
             npi.mods.all$npi_model[[2]],
             m_turnout_us13_all_npi,
             npi.mods.all$npi_model[[3]],
             m_turnout_us15_all_npi,
             m_midterm_us15_all_npi),
        custom.model.names = c("DK11 Participation",
                               "DK11 Turnout",
                               "US13 Participation",
                               "US13 Turnout",
                               "US15 Participation",
                               "US15 Turnout",
                               "US15 Midterm"),
        custom.coef.names = c("Intercept", "Narcissism (subset)",
                              "Female", "Age (2SD)", "Higher education", 
                              "Political interest (2SD)",
                              "Not Caucasian"),
        caption = "(Appendix) Narcissism, political participation, and turnout (with interest)",
        caption.above = TRUE)


screenreg(list(npi25.mods.all$npi_model[[1]],
               m_turnout25_dk11_all_npi,
               npi25.mods.all$npi_model[[2]],
               m_turnout25_us13_all_npi,
               npi25.mods.all$npi_model[[3]],
               m_turnout25_us15_all_npi,
               m_midterm25_us15_all_npi),
          custom.model.names = c("DK11 Participation",
                                 "DK11 Turnout",
                                 "US13 Participation",
                                 "US13 Turnout",
                                 "US15 Participation",
                                 "US15 Turnout",
                                 "US15 Midterm"),
          custom.coef.names = c("Intercept", "Narcissism (subset)",
                                "Female", "Age (2SD)", "Higher education", 
                                "Political interest (2SD)",
                                "Not Caucasian"),
          caption = "(Appendix) Narcissism, political participation, and turnout (25 NPI items, with interest)",
          caption.above = TRUE)

# 7-factor models, same order
screenreg(list(sub.mods$facet_model[[1]],
             m_turnout_dk11_facet,
             sub.mods$facet_model[[2]],
             m_turnout_us13_facet,
             sub.mods$facet_model[[3]],
             m_turnout_us15_facet,
             m_midterm_us15_facet),
        custom.model.names = c("DK11 Participation",
                               "DK11 Turnout",
                               "US13 Participation",
                               "US13 Turnout",
                               "US15 Participation",
                               "US15 Turnout",
                               "US15 Midterm"),
        custom.coef.names = c("Intercept", "Authority",
                              "Superiority", "Exploitativeness",
                              "Entitlement",
                              "Self-sufficiency",
                              "Vanity",
                              "Female", "Age (2SD)", "Higher education", 
                              "Exhibitionism",
                              "Not Caucasian"),
        caption = "(Appendix) 7-factors of Narcissism, political participation, and turnout",
        caption.above = TRUE)

screenreg(list(sub.mods.all$facet_model[[1]],
             m_turnout_dk11_all_facet,
             sub.mods.all$facet_model[[2]],
             m_turnout_us13_all_facet,
             sub.mods.all$facet_model[[3]],
             m_turnout_us15_all_facet,
             m_midterm_us15_all_facet),
        custom.model.names = c("DK11 Participation",
                               "DK11 Turnout",
                               "US13 Participation",
                               "US13 Turnout",
                               "US15 Participation",
                               "US15 Turnout",
                               "US15 Midterm"),
        custom.coef.names = c("Intercept", "Authority",
                              "Superiority", "Exploitativeness",
                              "Entitlement",
                              "Self-sufficiency",
                              "Vanity",
                              "Female", "Age (2SD)", "Higher education", 
                              "Political interest (2SD)",
                              "Exhibitionism",
                              "Not Caucasian"),
        caption = "(Appendix) 7-factors of Narcissism, political participation, and turnout (with interest)",
        caption.above = TRUE)

# 3-factor models
screenreg(list(sub25.mods$facet_model[[1]],
             m_turnout25_dk11_facet,
             sub25.mods$facet_model[[2]],
             m_turnout25_us13_facet,
             sub25.mods$facet_model[[3]],
             m_turnout25_us15_facet,
             m_midterm25_us15_facet),
        custom.model.names = c("DK11 Participation",
                               "DK11 Turnout",
                               "US13 Participation",
                               "US13 Turnout",
                               "US15 Participation",
                               "US15 Turnout",
                               "US15 Midterm"),
        custom.coef.names = c("Intercept", "Leadership/Authority",
                              "Entitlement/Exploitativeness",
                              "Grandiose exhibitionism",
                              "Female", "Age (2SD)", "Higher education", 
                              "Not Caucasian"),
        caption = "(Appendix) 3-factors of Narcissism, political participation, and turnout",
        caption.above = TRUE)

screenreg(list(sub25.mods.all$facet_model[[1]],
             m_turnout25_dk11_all_facet,
             sub25.mods.all$facet_model[[2]],
             m_turnout25_us13_all_facet,
             sub25.mods.all$facet_model[[3]],
             m_turnout25_us15_all_facet,
             m_midterm25_us15_all_facet),
        custom.model.names = c("DK11 Participation",
                               "DK11 Turnout",
                               "US13 Participation",
                               "US13 Turnout",
                               "US15 Participation",
                               "US15 Turnout",
                               "US15 Midterm"),
        custom.coef.names = c("Intercept", "Leadership/Authority",
                              "Entitlement/Exploitativeness",
                              "Grandiose exhibitionism",
                              "Female", "Age (2SD)", "Higher education", 
                              "Political interest (2SD)",
                              "Not Caucasian"),
        caption = "(Appendix) 3-factors of Narcissism, political participation, and turnout (with interest)",
        caption.above = TRUE)


# R version 3.5.1 (2018-07-02)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS  10.15.4
# 
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] texreg_1.36.23 polycor_0.7-9  psych_1.8.12   ggthemes_4.0.1 ggplot2_3.2.1  reshape2_1.4.3 broom_0.5.3    dplyr_0.8.3   
# 
# loaded via a namespace (and not attached):
#   [1] Rcpp_1.0.2       pillar_1.4.2     compiler_3.5.1   plyr_1.8.4       tools_3.5.1      zeallot_0.1.0    lifecycle_0.1.0 
# [8] tibble_2.1.3     nlme_3.1-137     gtable_0.3.0     lattice_0.20-35  pkgconfig_2.0.3  rlang_0.4.2      rstudioapi_0.10 
# [15] parallel_3.5.1   yaml_2.2.0       withr_2.1.2      stringr_1.4.0    generics_0.0.2   vctrs_0.2.1      grid_3.5.1      
# [22] tidyselect_0.2.5 glue_1.3.1       R6_2.4.0         foreign_0.8-70   tidyr_1.0.0      purrr_0.3.3      magrittr_1.5    
# [29] backports_1.1.5  scales_1.0.0     mnormt_1.5-5     assertthat_0.2.1 colorspace_1.4-1 stringi_1.4.3    lazyeval_0.2.2  
# [36] munsell_0.5.0    crayon_1.3.4  